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The behavior of the ground-state fidelity susceptibility in the vicinity of a quantum critical point 
is investigated. We derive scaling relations describing its singular behavior in the quantum critical 
regime. Unlike in previous studies, these relations are solely expressed in terms of conventional 
critical exponents. We also describe in detail a quantum Monte Carlo scheme that allows for the 
evaluation of the fidelity susceptibility for a large class of many-body systems and apply it in the 
study of the quantum phase transition for the transverse-field Ising model on the square lattice. 
Finite size analysis applied to the so obtained numerical results confirm the validity of our scaling 
relations. Furthermore, we analyze the properties of a closely related quantity, the ground-state 
energy's second derivative, that can be numerically evaluated in a particularly efficient way. The 
usefulness of both quantities as alternative indicators of quantum criticality is examined. 

PACS numbers: 75.10.Jm, 64.70.Tg, 03.67.-a, 02.70.Ss 
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I. INTRODUCTION 

The quantity known as fidelity naturally appears in 
the field of quantum information science as a way of de- 
termining the reliability of a given protocol for quan- 
tum information transfer: the similarity between input 
| 'I'in) and output |* ut) states can be quantified by sim- 
ply computing the absolute value of the overlap be- 
tween them, F = |(\&i n |^out)|- Recently, after the pi- 
oneering work 1 - of Zanardi and Paunkovic, and following 
the broader trend of cross-fertilization between the fields 
of quantum information science and condensed matter 
physics^ a number of studies have extended the scope 
of applicability of the concept of fidelity to the study of 
quantum critical phenomena (for a review, see Ref. |3|). 

The basic idea behind this so-called fidelity approach is 
simple. We consider a general many-body Hamiltonian 



H(g) = H +gH 1 



(1) 



with ground-state |* (flO>, ft(ff)|#o(ff)) = E (g)\^ Q (g)). 
Since | v f'o(5)) undergoes major changes in the vicinity of 
a quantum critical point (QCP) g c , we expect a sharp 
drop in the fidelity, 



F(g,dg) = \{* (g + dg)\* (g))\ , 



(2) 



for small (dg — > 0) variations in g close to g c . There- 
fore, by investigating the behavior of F(g, dg) when cou- 
plings in the Hamiltonian are varied, one should be able 
to detect quantum criticality. Besides its novelty, this 
approach is purely quantum geometrical^ and therefore 
has the appeal that no a priori identification of order 
parameters is required. 

The concept of fidelity susceptibility^ Xf(<?) naturally 
appears as the fidelity's leading term in the limit dg — > 0, 



F(g,dg^0) 



1 - ]pCF{g)dg 2 . 



(The linear term in dg in the above expansion vanishes 
due to normalization of the wave-function — alterna- 
tively it can be seen to arise from the fact that F(g,dg) 
is maximum at dg = for any value of <?.) The aforemen- 
tioned drop in F(g, dg) close to a QCP is thus associated 
to a divergence in xf and the latter quantity may also 
be employed in the study of quantum phase transitions. 
The situation here is reminiscent of the use of the spe- 
cific heat to detect thermal phase transitions: while the 
presence of singularities in the specific heat for varying 
temperatures signals the location of finite-temperature 
critical points, xf(<?) is a system's response to changes 
in the coupling constant g, whose divergencies are asso- 
ciated to the occurrence of quantum phase transitions. 

Although obviously some information is lost in going 
from F(g, dg) to Xf(<?), an d for instance it is currently not 
clear whether transitions of order higher than second can 
be detected by studying the latter, focusing on \p(g) lias 
up to now proved to be a fruitful strategy. The main rea- 
son behind this is that it is possible to show£~— that Xf(<?) 
is closely related to more conventional physical quanti- 
ties, such as imaginary-time dynamical responses. This 
is particularly advantageous since it allows one to rely on 
well established concepts and techniques from theoretical 
condensed-matter physics in order to draw conclusions on 
the properties of xf(<?)- We follow this line of reasoning 
in this paper, in a twofold way. 

First, we present the details of a recently introduced^ 
quantum Monte Carlo (QMC) scheme that allows for the 
evaluation of xf for a large class of sign-problcm-frcc 
models. This constitutes an important advance as the 
group of problems that can be studied within the fidelity 
approach is considerably enlarged, and additionally one 
benefits from the computational power of QMC methods. 
In particular, high-precision scaling analysis for models 
in dimensions higher than one is now possible: previous 
computations of xf for two-dimensional systems have re- 
lied on exact diagonalization (ED) techniques and were 
restricted to small system sizes, something that precludes 
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a precise determination of scaling dimensions in the vicin- 
ity of a QCP. 

Second, by building upon the aforementioned relation- 
ship between xf and response functions, we determine 
the scaling behavior of the fidelity susceptibility close to 
a QCP. The divergence of xf(<?) at g c is shown to be re- 
lated to the critical exponent v describing the divergence 
of the correlation length. In this way, and supported by 
the results obtained from the QMC simulations, we assess 
the validity of other scaling analysis for the divergence of 
Xf that have recently appeared in the litterature. 

Throughout the paper, we also analyze the prop- 
erties of the ground-state energy's second derivative, 
d 2 Eo(g)/dg 2 . This quantity is closely related to the fi- 
delity susceptibility^ and, as explained in Sec. Ill Cl scales 
in a way related to the scaling of xf in the quantum criti- 
cal regime. Since the computation of d 2 Eo(g)/dg 2 within 
QMC is much more efficient than that of xf, as discussed 
in Sec. MI B[ it is important to address the question of 
which of these quantities is best suited to the study of 
second-order quantum phase transitions and of whether 
the current interest around the concept of fidelity suscep- 
tibility is justified on practical grounds. 

The paper is organized as follows. After reviewing ba- 
sic concepts in the fidelity approach in Sec. Ill A\ we ana- 
lyze an extension of the concept of fidelity susceptibility 
to finite temperatures (a prerequisite for path-integral 
QMC simulations) and relate it to a more commonly em- 
ployed metric for thermal states in Sec. Ill Bl We then 
perform a scaling analysis of xf and d 2 Eo(g)/dg 2 in 
Sec. Ill CI and relate their scaling dimensions to conven- 
tional critical exponents. 

In Sec. IIIII we give a detailed account of the previ- 
ously introduced^ QMC scheme for calculating xf and 
further explore it (Sec. IIV[) in order to determine the 
scaling dimension of the fidelity susceptibility in one of 
the most paradigmatic models in the field of quantum 
phase transitions: the transverse-field Ising model (TIM) 
in two dimensions. Throughout the paper, concepts are 
illustrated by presenting results for the one-dimensional 
version of the TIM, which is exactly solvable^— A sum- 
mary is given in Sec. |V] and important technical de- 
tails are discussed in the Appendix. Some of the results 
present in thispaper have been first presented (by some 
of us) in Ref. 



II. FIDELITY SUSCEPTIBILITY 



Definition 



We consider the limit of dg — > and perturbatively 
calculate the overlap appearing in Eq. ^ to leading order 
in dg. The fidelity susceptibility, defined by Eq. (JTJ) , can 
easily be shown to read^^ 



Xf(sO 



E 



[E n (g) - E (g)} 2 



(3) 



in terms of the eigenbasis 

£|*„(s)><*»(ff)l = -r (4) 

n 

oiH(g), H(g)\* n (g)) = E n (g)\* n (g)). 

Starting from Eq. ([3]), one can relate^ xf(s) to the 
imaginary-time correlation function 



G Hi (t)=6(t)((H 1 (t)H 1 (0))-(Hi) 2 ) 



(5) 



where H\(t) = e T ' H Hie~ T ' H , with r denoting an 
imaginary time, (9(r) is the Heavisidc step function 
and zero-temperature averages are defined by (O) = 
(^ (g)\O\y (g)) . Inserting Eq. Q into this last equa- 
tion and taking its Fourier transform we arrive to 

- \ - |(vMg)|gil*o(g))l 2 

^ E n (g) - E {g) + iuj 



The similarity between Eqs. ^ and ([5]) is evident and by 
simply performing a derivative, one can establish^ the 
important result 



Xf(s) 



dto 



j=0 



dTTG Hl (r) . (7) 



This expression is remarkable for a number of reasons. 
First, it relates xf(<?) to a dynamical response of the 
system to the "driving term" Hi, evidencing its physi- 
cal content. Second, as discussed in detail in Sec. Ill CI 
Eq. ([7]) allows us to address the issue of the scaling be- 
havior of Xf(<?)- Finally, Eq. ([7]) permits us to extend 
the definition of fidelity susceptibility to finite temper- 
atures (Sec. Ill B[) and therefore constitutes an obvious 
starting point in devising a scheme to obtain xf(s) from 
path-integral QMC simulations. 

Before proceeding, it is also instructive to consider the 
ground-state energy's second-derivative, whose intimate 
relation to xf(<?) has been pointed out in Refs. p7l andfT^. 
Motivated by this close relationship, we define xe(<?) = 
—d 2 E (g)/dg 2 and, by using the eigenbasis of H(g), it is 
readily shown that 



Xe(s) 



d 2 E (g) 
dg 2 



E n (g)-E (g) 



(8) 



Comparing Eqs. (j6|) and (jSJ) it is straightforward to see 
that 

/>oo 

Xe(9) = 2G Hl (w = 0) = 2 / drG Hl (r) . (9) 

Jo 

We can notice that the only important difference between 
Eqs. © and ([5]) is that in the former the denominator is 
squared: this is reflected by the appearance of the r fac- 
tor in Eq. (0, absent in Eq. ©. One might thus expect 
Xf(<?) to display a more pronounced behavior around a 
QCP and therefore to be a better indicator of quantum 
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criticality, an observation put onto firmer grounds by the 
scaling analysis of Sec. Ill Cl Finally, similarly to the case 
of Eq. ([7]), the relation Eq. (j9|) can be used in order to 
extract Xe(sO from QMC simulations, as we discuss in 

sec. hub 



B. Finite Temperature 

Before discussing how to extend the definition of fi- 
delity susceptibility xf(<?) to finite temperatures (T = 
1//3) it is instructive to consider first the similar exten- 
sion for xe(s)- From Eq. (|9]), one obtains the finite-T 
generalization 

/■/V2 

Xe(<7,/3)=2/ G Hl (r)dT , (10) 
Jo 

where now Gh-i (r) is still defined by Eq. ([5]) but with 
thermal averages, (O) = Z _1 Tr [cxp(— f3W)0] replacing 
ground-state expectation values [Z = Tr {exp(— (3H)} is 
the partition function] . An important subtlety is appar- 
ent here: notice that the upper integration limit in the 
above expression is /3/2, instead of (3. The underlying 
reason is that, within the path-integral formalism used 
in QMC simulations, periodic boundary conditions are 
implied along the imaginary time direction. Connected 
physical correlation functions, such as Ghi{t), are pe- 
riodic along the T-direction, with period f3. This is il- 
lustrated in Fig. [IJa) for the one-dimensional TIM (see 
Sec. lIV A| on a chain with L = 16 sites and (3 J = 16, for 
h/J = 1: we see that Gh x {t) is symmetric around (3/2 
(vertical dashed line) and decays to zero at r —¥ (3/2 for 
large enough f3, a trend already noticeable in Fig. Ufa) 
where data for the relatively high temperature (3 J = 16 
are displayed. Therefore, in the present case we have 
XE (g) = 2 fj 2 G Hl (r)dr = G Hl (r)dr. 

The definition of fidelity susceptibility [Eq. ([7])] can be 
extended to finite temperatures in a similar way 

Xf(<?,/3)= / tG Hi {t)cLt . (11) 

JO 

An important difference appears though: the aforemen- 
tioned properties of Gn t (t) are not shared by the func- 
tion tGrx (t) , since the pre- factor r destroys the pe- 
riodicity along the imaginary-time direction This is il- 
lustrated in Fig. [T](b), again using the one-dimensional 
TIM as an example. In particular, J Q tGhx (t)<2t 7^ 
Ip/2 t Ghx {t)cLt and, therefore, in order to ensure that 
Xf(<7, (3) converges correctly to its zero temperature limit, 
lim^oo Xf{9, (3) = XF(ff); one must cut the integral at 
(3/2. This has important implications for the QMC eval- 
uation that is now made possible by Eq. (fTTj) , as clarified 
in Sec. HITAl 

While the just discussed generalization of \f to finite 
(3 [Eq. (|11[) ] has been introduced for computational pur- 
poses, it is possible to relate it to more commonly used 




t J rJ 

FIG. 1: (Color online) Correlation functions (a) Gh 1 {t) and 
(b) tGh x (i~), as a function of r, for finite inverse temperature 
/3J — 16 for the d = 1 TIM on a chain with L = 16 sites 
and h/J — 1 (see Sec. IIV A[) . Data are generated by using 
the QMC method detailed in Sec. lIIII and fixing the parity to 
the P = +1 sector (Sec. lIV A~2|) . Exact Diagonalization (ED) 
data with fixed parity P = +1 are also shown both for finite 
and zero temperatures. 

metrics for thermal (mixed) quantum states^ as we dis- 
cuss in what follows. 



1. Bures Metric 

The so-called Uhlmann fidelity generalizes the concept 
of fidelity [Eq. @] to the case of mixed states. For den- 
sity matrices /?a and pb, it is defined as^ 



F(pa,Pb) = Try p^PBp 1 / 2 . 

and has an associated metric ds(pA,PB) = 
\J2 [1 — F(pa, Pb)], known as the Bures distance. 
We arc interested in the case of thermal density 
matrices, 

p g = ^J2 e ^ Enig) \^mn(g)\ , 

n 

expressed here in terms of the eigenbasis oiW(g) [Eq. (j4])] . 
The concept of fidelity susceptibility [Eq. ((3])] can then 
be extended to the finite-T regime with the Bures metric 
ds 2 (g, /3) = ds 2 (p g , p g +d g ) for density matrices associated 
to infinitesimally close {dg — > 0) couplings in the Hamil- 
tonian, g and g + dg. Following Zanardi et al.^ one ob- 
tains the following expression for ds 2 (g,{3) (for the sake 
of simplicity we omit the dependence on g of the eigen- 
values and eigenvectors in the remainder of this Section) 

ds 2 ( 5 ,/3)=d4( ff ,/3) 

V [(^[^[^[^-^(l-e- 2 -) 2 (12) 
{E n -E m f Z l + e- 2 * ' 
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where x = f3(E m — E n )/2. One appealing feature in this 
expression is that it distinguishes "classical" and quan- 
tum contributions. The "classical" term, ds 2 x (g,f3), is 
given by^ 

d S 2 cl (g,f3) = ^((Hl d )-(H hd r) . (13) 

Here, Hid denotes the diagonal elements of H\ in the 
cigenbasis Eq. (jH), so that 



On the other hand, the second term in Eq. (|T2j) is of pure 
quantum origin and vanishes unless [p g , pg+dg] 7^ 0.— 



2. Relation between xe{g, ft) and ds 2 (g , /3) 

We now relate the two previously discussed finite tem- 
perature extensions for the fidelity susceptibility, namely 
Eqs. (|TT|) and ([THUS])- We expand the trace for thermal 
averages and insert the eigenbasis of T-L(g) [Eq. flU] in 
Eq. (fTTj) arriving to 



Xf(9>P) 



Z^ 



0/2 
d,TT 



-PE n +T(E n -E„ 



K*„|ffl|*, 



(E„+E m ) 



(*„|ffi|# n )<* m |ffi|* r 



The terms with n = m in the first term in the integrand 
are independent of r and can be regrouped with the sec- 
ond term to yield ^ds^(g, /3). Performing the integration 
for the remaining terms, we finally obtain 

v (a ff\ - da cl(g'^) 

Xf{9,P) - 



E 

71 > 771 



I 2 e-f* 15 " 



(14) 



(E n — E„ 



where again we set x = (3(E m — E n )/2. One can readily 
show that in the limit T -> both Xf(9,P) [Eq. ([Ti])] 
and ds 2 (g,f3) [Eq. (JT2J)] converge to the ground-state re- 
sult of Eq. ([3]), as desired. This is illustrated for the 
TIM on the square lattice in Fig. [3l where data from 
QMC [for xf (<?,/?)] and exact diagonalizations [for both 
Xf(Sj/3) and ds 2 (g,/3)] are displayed (see discussion in 
Sec. HVBp . On the other hand, the high-temperature 
limit (ft -> 0) yields ds 2 (g,f3 -> 0) = 2 X f{9,P -> 0) = 
f ((H 2 ) - ( Hl ) 2 ). 



C. Scaling Behavior 

We focus now on the issue of how the fidelity suscep- 
tibility \F behaves in the vicinity of a QCP and start 
by reviewing the main results from the scaling analysis 
devised by Campos Venuti and Zanardi.— 

Starting from Eq. ([7]) and following standard scaling 
arguments (see for instance Ref. [H) , they apply the scale 
transformation x' = sx, t' — s z t (z is the dynamic criti- 
cal exponent), and arrive to the following relation for the 
scaling of the fidelity susceptibility density: 



Xf ~ \g 



ii/(d+2z-2A ffl ) 



(16) 



In order to analyze the general case, we evaluate the Hi 



Here, L d = N is the number of sites for a d-dimensional 
system. In what follows, we assume that a second-order 
quantum phase transition takes place at a value g c of the 
"driving parameter" g and that the correlation length 
diverges in its neighborhood as £ ~ \g — g c \~"\ is 
the scaling dimension of the driving term H\ in Eq. |T]): 



ratio between f(x) = (1— e 



)7(i+ 

in Eq. dHJ), and g(x) = (1 - e~ x f [from Eq. ([14])]. We 
have f(x)/g(x) = 1 + 1/ cosh(a;) and therefore f(x)/2 < 
g(x) < f (x). Noting that ds 2 cl {g,P) > 0, we can conclude 
that 

\ds 2 {g,(i) < X f(9,P) <ds 2 (g^) . (15) 

These inequalities show that if ds 2 (g,ft) diverges, then 
X p(g,P) must also diverge and we conclude that both 
quantities are equally well suited for detecting criticality. 
While in this paper we are interested in quantum phase 
transitions and thus focus on the limit (3 — > 00, it would 
also be of interest to investigate the finite-T behavior of 
Xf{9iP)'- as discussed in Ref. [H, this quantity might be 
able to detect thermal phase transitions and/or finite-T 
signatures of quantum criticality. 



"), that appears thus imply that 



H ^H\. Standard finite-size scaling arguments 



XF 



-2(z-A Hl ) 



(17) 



for finite systems at criticality. 

An alternative scaling analysis has been recently put 
forward by some of us in Ref. [8|, where the simpler result 



XF 



L 2/u 



(18) 



has been derived. The appeal of the above scaling rela- 
tion, as compared to Eq. (fTT)) . stems from the fact that 
it is not expressed in terms of the exponent , rarely 
dealt with in more conventional approaches to quantum 
critical phenomena. In fact, A/^ can easily be related to 
the exponent v (see Ref HI). In what follows, we provide 
an intuitive equivalent derivation of the result Eq. (|18[) . 
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We start by remarking that at T = the free-energy 
density /h reduces to the ground-state energy density, 
L~ d E . Since, by definition, /h ~ \g ~ ffc| 2_a (see for 
instance Ref . [la ) , one readily obtains 



T -d 



Is - .9c 



(19) 



This relation is not at all surprising since, as we can 
see from Eq. ©, Xe(<?) = —d 2 E {g)/dg 2 is similar to a 
"zero-tcmpcrature specific heat" . Comparing the expres- 
sions for xe [Eq. ©] and xf [Eq. (0)], one sees that the 
only difference is the presence of an imaginary time (or 
inverse energy) scale r in the latter, that we naturally 
expect to scale as \g — g c \~ zv m the critical regime^ We 
thus arrive at 



XF 



9c 



-{a+zv) 



= \g 



-(2-dv) 



(20) 



where we have made use of the hyper-scaling formula 
2 — a = v{d + z). Therefore, we arrive to the conclusion 
that L~ d XF diverges only when v < 2/d. 

On the other hand, by similarly inserting the hyper- 
scaling formula into Eq. ([TO)) , we obtain L~ d \E ~ Iff — 
g c \- 2 +v( d + z ) and, as anticipated in Sec. Ill A[ conclude 
that xe has a weaker divergence than xf at a critical 
point. Furthermore, wc find that the more stringent con- 
dition v < 2/(d+ z) must be satisfied for L~ d \E to 
diverge and that there might be situations where only 
L~ d \F displays a divergence at a critical point, being in 
general a better indicator of quantum criticality. 

Next, by performing a finite-size scaling analysis, we 
conclude that for the fidelity susceptibility per site we 
have 



XF 



L> 



in agreement with the result Eq. 
L~ d \E, wc similarly obtain 

L~ d XE ~ Li- {d+z ^ 



(21) 

W) (Ref. i). For 



(22) 



function^ 2 - 



Z = Tr (e~ m ) ~ 



E E E ^7 

n=0 a S„ 



(23) 

Here {\a}} is any suitable basis and the system's Hamil- 
tonian is typically a sum over local operators: H = 
'}2n,T~L{b), with b labeling different local terms. For in- 
stance, b may denote operators acting on different bonds 
of the lattice and/or diagonal versus non-diagonal opera- 
tors. For our current purposes, it is convenient to choose 
a decomposition that respects the bipartition of Eq. (fTJ) , 
such that all terms appearing in Hq are labelled by bo and 
those appearing in gHi by b\ and we have b € {bo,bi}. 
SSE configurations (a,5„), with operator strings 



S n = f[H(b i ) , 



(24) 



are then sampled, according to the statistical weight 



W{a,S n ) 



; = 1 



Efficient update schemes such as the directed loop 
algorith m 21 ' 23 render the SSE technique one of the most 
efficient QMC methods for quantum lattice models. 

The general procedure for obtaining thermal averages 
within the SSE framework is discussed in detail by Sand- 
vik in Ref. |24| . The basic idea, supposing we are inter- 
ested in an observable 0, is to determine an estimator 
0(a, S n ) such that 



E 0(a,S n )W(a,S n ) . 



« (a,S„) 



In what follows, we show how estimators for the fidelity 
susceptibility XfO,£) [Eq. dTTJ] and Xe(.9,/3) [Eq. (fTU]>] 
can be obtained from SSE QMC simulations. 



The validity of the scaling relations Eqs. (|2Tj) and ([221) is 
confirmed by the analysis of the QMC data performed in 
Sec. lIVCl Finally, we note that the scaling relations pre- 
sented here have been independently found in the context 
of quantum quenches ^r— as discussed in Sec. [Vj 



III. STOCHASTIC SERIES EXPANSION 

Despite the sign problem that precludes efficient sim- 
ulations of most fermionic/frustrated models, Quantum 
Monte Carlo methods are among the most efficient tools 
to simulate quantum many-body problems. Particu- 
larly useful for our purpose here is the QMC formula- 
tion known as Stochastic Series Expansion (SSE), devel- 
oped by Sandvik and coworkers*^— Basically (for a de- 
tailed account the reader is referred to Ref. [211 ) , SSE re- 
lies on a power series expansion of the system's partition 



A. Fidelity Susceptibility 

First, we need to evaluate imaginary-time operator 
products of the form (Hi(t)Hi(Q)) appearing in the inte- 
grand of Eq. pT[) [cf. Eq. (J5J]. These operators being part 
of the Hamiltonian, one trick consists in re-interpretating 
two of the elements with label b\ of the string E q. (12"4"1) 
as the operators to be measured. Following Ref. [24|, we 
arrive to 

g 2 (Fi(r)F!(0)) = 

n-2 



E 



(n-l) 



— ' (n — m — 2)!m! 



T) n-m-2 T m ( NgHl {m)) w . 



(25) 



Here, n is the length of the operator string S n [Eq. 
and N g H 1 (m) the number of times any two operators 
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comprising gli\ appear in the strings S n separated by 
m positions. We discuss below how N g jj 1 (jn) can be 
measured. 

The second term in Eq. (JS|) is obtained by a simpler 
procedure^ and is given by 



(Hi) 



1 



g 2 P 



(N gHl : 



w 



(26) 



where N g H x is the total number of gli\ operators in S n . 

Inserting the results Eqs. (|25l l26|) into Eq. (fTTj) and 
integrating from r = to j3/2 (taking into account the 
important multiplicative factor of r in the integrand), we 
finally arrive to the result 



n—2 



I " " 

Xf(9,P) = ~2 XI [ A (™, n ) (%iH)wl 



g 

a m=0 

with the coefficient 



(27) 



A(m, n) 



(n-l)l 



(n — m — 2)!m! 



1/2 



-Hi-rr 



(28) 

We show in Appendix [S] how this coefficient can be ap- 
proximated very accurately by an analytical expression 
in the limit of n^$> 1. 

NgHx ( m ) is conveniently extracted from the simula- 
tions in two steps. Firstly, the string Eq. (j24|) is tra- 
versed (for instance when performing diagonal updates; 
see Ref.Qjl) and the positions i where a local Hamiltonian 
H(b l ) appears with a label b l — b\ are recorded (there are 
in total N g H 1 such operators). Secondly, the histogram 
NgH-y (tti) is generated by computing all distances m be- 
tween all previously recorded positions i. This step is the 
most demanding as it requires N g H 1 {N g H 1 — l)/2 oper- 
ations. Note finally that the prefactor 1/g 2 arises from 
the definition of the fidelity susceptibility Eq. ([3]) which 
docs not include the coupling constant g, whereas the 
SSE decomposition used in Eq. (|2"5|) typically does. 



B. Ground-State Energy's Second Derivative 



The results Eqs. (|251 126| can also be used in order to 
directly evaluate the ground-state energy's second deriva- 
tive, relying on Eq. (|10[) and extrapolating to the limit 
f3 — > oo. The absence of the factor t in Eq. © consid- 
erably simplifies the situation since the integration over 
t can now always be performed exactly. In this way, we 
arrive to the simple result 
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XeG?,/?) = [(N 2 gHl ) w - (N gHl ) w - (JW; 

(29) 

We stress that the computational cost for evaluating 
Xe (<?,/?) is much lower than the one required to ob- 
tain XF{g,fi)'- the estimator for the former quantity in 
Eq. (|29[) simply requires counting the number of times 
the operators contained in the "driving term" gli\ occur 



in the operator strings S n . This is to be contrasted with 
the computationally heavy task, specially in the limit of 
large lattice sizes and low temperatures, of computing 
the histogram NgH^rn) necessary in evaluating Xf(<?>/3) 
[cf. Eqs. ma ESM- 



IV. NUMERICAL SIMULATIONS 
A. Transverse-Field Ising Model 

1. Definition 

The transverse-field Ising model (TIM) is perhaps the 
simplest model to display a QCP and many key concepts 
in the theory of quantum critical phenomena have been 
developed by analyzing its properties^ 

The TIM Hamiltonian reads 



H(h) = JHj + hH h 



i 3 



(30) 



where denotes nearest-neighbor sites on a d- 

dimensional lattice and a^' z are Pauli matrices attached 
to the site i. We set the energy scale by henceforth fixing 
J = 1. 

At zero temperature and in any dimension, a quantum 
phase transition occuring at a field h c separates a ferro- 
magnetic phase for low fields h < h c from a polarized one 
for h > h c , with spins aligning along the field direction. 
The quantum phase transition in d dimensions belongs 
to the universality class of the finite-temperature phase 
transition of the classical Ising model in d+ 1 dimensions 
and has a dynamic critical exponent z = l3& 

The TIM is exactly solvable in the one-dimensional 
case^— The QCP is located at h c = 1 and most observ- 
ables, including XFr can be computed analytically. We 
make use of these exact results in establishing the valid- 
ity of the QMC method discussed in Sec. lIIIl for instance 
m Sec. IIVA21 where the issue of parity is discussed. On 
the other hand, the TIM is not solvable in two dimen- 
sions and has been investigated mainly through means 
of numerical techniques^ - — The most precise estimate 
for the location of the QCP, h c = 3.04438(2), has been 
obtained from a QMC approach^ 

Before proceeding, we remark that a scaling analysis 
for the fidelity susceptibility and the second derivative of 
the ground-state energy for the TIM on the square lattice 
has been recently performed by Yu and collaborators in 
Ref. |H. Results are compared in Sec. El 



2. Parity quantum number 

An important issue concerning the TIM is the exis- 
tence of a conserved quantum number, the parity P. It 
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FIG. 2: (Color online) Fidelity susceptibility density for the 
Id TIM, as obtained from the exact solution (curves) and 
QMC simulations (symbols; see Sec. IIII A|) . Results for both 
parity sectors P = — 1 and +1 are shown (see Sec. IIV A 21) . 
QMC results indicated by P = ±1 (triangles) have been ob- 
tained by making no distinction between parity sectors. For 
h/J > 1, statistics are insufficient to estimate xf m the 
P = — 1 sector from QMC simulations since only the ground- 
state in the P = +1 sector is sampled at the low temperature 
considered here, /3 J — 32. Data are for a system of size L = 16 
with periodic boundary conditions. 

is readily verified that the parity operator, 

N 

^ = ll a i ( 31 ) 

i=l 

(for a system with N sites), commutes with the TIM 
Hamiltonian Eq. P"U)) , [H,V] = 0, and therefore the par- 
ity P = ±1 is a good quantum number. 

For finite systems, the ground-state of the TIM lies in 
the P = +1 sector, as shown by the following argument. 
It is convenient to work in the basis given by tensor prod- 
ucts of the eigenvectors {11%, |-J%} of of at every site: 
{\(f>m}} with, for instance, \</> m ) = |tWt ••%• All off- 
diagonal matrix elements for the Hamiltonian Eq. (|30[) 
are non-positive in this basis and therefore the TIM 
on finite lattices satisfies the conditions for the Perron- 
Frobenius theorem to apply. According to this theorem, 
the coefficients of the system's ground-state (in its expan- 
sion in terms of the basis {!</>,„)}, |*o) = XL C ™I^™)) 
must all have the same sign (say, c m > 0). Consider the 
lowest-lying states in each parity sector \^q)- they can 
be expanded as |<Fj) = £ m c±|0±), with {1^™)} denot- 
ing the subset of elements in {|^ m )} with fixed parity 
P = ±1. The argument proceeds by noticing that the 
parity operator [Eq. ([31]) ] simply acts as a spin rever- 
sal operator upon the elements of the a x basis: namely, 
^l^m) = IV'm). where l^m) is obtained from |</>±) by flip- 
ping all spins (for instance, P|t44-t • • % = \lftl ■ ■ %)• 
Thus, the expression 

m m 



is only consistent with the positiveness of the ground- 
state if P = +1 (one readily sees that the above relation 
can only be satisfied in the P = — 1 sector if the coeffi- 
cients for the basis elements and V\4>„) = \ipn) have 
opposite signs in the expansion for )). 

The above discussion is directly relevant for our pur- 
poses here since, while expectation values for most phys- 
ical observables are the same for the lowest-lying states 
in both parity sectors, this turns out not to be the case 
for xf and %e- This is illustrated in Fig. [21 where exact 
results for xf (curves, see Ref. 0) for the TIM on a chain 
with L = 16 sites are shown for both parities, P = +1 
and P = — 1. Also shown are data obtained from a naive 
QMC implementation not discriminating between differ- 
ent parity sectors (triangles in Fig. [5]): the disagreement 
is evident, specially in the neighborhood of the QCP at 
h c = 1. While we expect this discrepancy to disappear 
in the thermodynamic limit (as suggested by exact re- 
sults for the d = 1 TIM for increasingly larger systems; 
however, we have no general proof), QMC simulations 
are obviously restricted to finite system sizes and it thus 
important to take the parity quantum number into ac- 
count. 

Fortunately, the parity of a given SSE configuration 
(a, S„) can easily be determined within the here adopted 
convention for the TIM [Eq. ([50])]i2£ Indeed, the parity 
operator, defined in Eq. pip , is diagonal in the a z basis 
employed in SSE-QMC simulations and the parity is then 
readily obtained as 

AT 
i=l 

(since the parity is a conserved number, it can be com- 
puted at any time-slice, in particular at r = 0). Both 
parity sectors are sampled due to the non-local updates 
in the QMC scheme (see Ref. |2T[ ) . Therefore, our strategy 
consists in computing the estimators required for obtain- 
ing xf and Xe (see Sec. IIII|) for all SSE configurations 
(a, S n ) and storing the results in different variables ac- 
cording to the parity of the state \a(r = 0)). Data ob- 
tained in this way for the d = 1 TIM are shown in Fig. [5] 
and perfectly agree with the exact results. Finally, since 
for finite systems the ground-state has P = +1, all re- 
sults discussed in the present work have been obtained 
for this parity sector (with the exception of the indicated 
ones in Fig. [2]). 

B. Simulation Details 

The computation of xf (<?,/?) and xe(<7,/3) requires 
only small changes to an existing SSE code: estima- 
tors for both quantities [Eqs. (|27|) and ([%!?]) ■ respectively] 
are simply computed by analyzing the operator strings 
[Eq. (|24|) ]. a task ideally carried out while performing di- 
agonal updates for the SSE configurations (a, S n )^- Our 
code is based on the ALPS^i librairies implementation 
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FIG. 3: (Color online) Finite-T fidelity susceptibility 
L ~ 2 XF{g,/3) [Eqs. H11I14|> ] and Bures metric L~ 2 ds 2 (g, /3) 
[Eq. (| 12f) ] as a function of the inverse temperature /3 = T -1 for 
the TIM on the square lattice (iV = 4x4 sites cluster with pe- 
riodic boundary conditions; restriction to the P = +1 parity 
sector, see Sec. lIV A2l) . Both Hj and hH h [see Eq. ([30])] have 
been considered as the "driving term" Hi in the definitions for 
Xf (<?,/?) and ds 2 (g,/3): we see that different choices yield the 
same T = result. The vertical dashed line marks the value 
j3 = 2L set in obtaining the QMC data displayed in Figs. [2] 
and 21 while for the coupling considered here Ql/J = 1) the 
system is deep into the ferromagnetic phase and convergence 
to ground-state expectation values is achieved for smaller /3, 
the more stringent condition f3 — 2L is necessary closer to the 
QCP. 



of SSE QMC^ The main modifications of the original 
codes are independent from measurements for xf and 
Xe and are specific to the TIM [as defined by Eq. (f3T)]) ] 
studied in the present work. They involve changes in the 
processes that are allowed when performing off-diagonal 
updates^ and the computation of the parity quantum 
number (see Sec. IIV A~2j) . 

Calculating xf can be computationally demanding due 
to the fact that the number of operations required for 
obtaining the histogram NgH^m) [sec Eq. ([27])] scales 
quadratically with the total number of gH\ operators 
in the string S n . The situation can be ameliorated by 
a judicious bipartition of the system's Hamiltonian into 
Hq and H x [Eq. ©]. 

Indeed, there is freedom to consider either Hj or hHh 
(or, in general, any combination of these) appearing in 
Eq. ((301) as the "driving term" H x in Eq. fl5J): since 
(V n {h)\(Hj + hH h )\V (h)) = for n ^ 0, we readily 
conclude from Eq. ([3]) that Xf J = h 2 XF Hh (superscripts 
indicate the term assigned to Hi) and therefore different 
bipartition choices lead to the same zero-temperature re- 
sults for xf, apart from a trivial multiplicative factor. 
Although this is no longer true for the finite-T general- 
izations of Eqs. (fl"2j) and (|14p .— the equivalence between 
results obtained from different partitions is recovered in 
the limit T — > 0, as shown in Fig. [3] for the TIM on the 
square lattice (data obtained from QMC simulations for 



a 4 x 4 cluster and h/J = 1). 

We may thus explore the fact that different terms 
in the Hamiltonian dominate in different regions of the 
phase diagram in order to reduce computational cost. 
Specifically, for the case of the TIM considered here 
[Eq. ([30])] it is more efficient to compute N g jj 1 (m) in the 
high- field limit if we set H\ — Hj, since the number of 
such operators in the strings S n will be smaller than that 
of Hh operators in this limit. In practice, we find that, 
for field magnitudes close to the QCP for the TIM on 
the square lattice, setting Hi = Hj is the most efficient 
choice. 

Following this strategy, we have simulated the TIM on 
the square lattice by considering clusters with linear size 
L and periodic boundary conditions (PBC), and are able 
to reach L = 28 when computing xf- On the other hand, 
computing xe requires much lesser numerical effort and 
we are able to reach L = 48. For both quantities, we find 
that if we set the inverse temperature (3 = 2L both xf 
and Xe reach their ground-state expectation values. This 
is illustrated in Fig. [3] for the L = 4 cluster. We have also 
performed a few simulations setting /3 = 4L in order to 
confirm that convergence has indeed been achieved, at 
least within error bars, for /3 = 2L. 



C. Results 

Our QMC data for L~ 2 xf and £~ 2 xe for the TIM on 
the square lattice are shown in Fig. 2] for various system 
sizes L. The presence of peaks in the curves for both 
quantities is evident: they become more pronounced for 
increasing L and their positions seemingly converge to- 
ward the estimate h c — 3.04438(2) for the QCP found in 
Ref. l3ll (see below). Furthermore, we notice that xe dis- 
plays less pronounced peaks than xf, as expected from 
our discussion in Sec. Ill CI A quantitative data analysis 
is explained in what follows. 

We start by determining the peaks positions and 
heights for both xf and xe from the raw data displayed in 
Fig.fU The so obtained results are shown in Fig. [S] From 
the scaling relations derived in Sec. Ill Cl L~ d XF ~ L»~ d 
and L- d XE - [Eqs. ([2l]) and ([22]): d = 2 and 

2 = 1], we expect a linear dependence for the logarithm of 
the peaks' height on In L. This is confirmed by the results 
shown in Figs. [SJa) and (c). By applying linear regres- 
sion to the points associated to the three largest values 
of L in each plot we obtain our first estimates for cor- 
relation length's exponent: v = 0.623(8) [xf, Fig. EJa)] 
and v = 0.615(1) [xe, Fig. [SJc)]. While the former esti- 
mate is in good agreement with the result for the univer- 
sality class of the three-dimensional classical Ising model 
[v = 0.6301(8), Ref. [I?}, the latter clearly underestimates 
v. This is likely to be explained by the weak divergence 
displayed by xe, implying that regular sub- leading cor- 
rections are important in accounting for the behavior in 
system sizes as the ones considered here: indeed we no- 
tice that the data points corresponding to the smallest 
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FIG. 4: (Color online) (a) Fidelity susceptibility density 
L~ 2 \¥ and (b) ground-state energy's second derivative per 
site L~ 2 d 2 E (g)/dg 2 = -L~ 2 xv{g) for the TIM on the square 
lattice, as a function of h/J and for indicated system sizes L 
(temperatures are set to f3 = 2L). Data have been obtained 
by applying the SSE QMC procedure detailed in Sec. IIIII 



system sizes clearly deviate from the linear fit obtained 
for the points for the three largest L in Fig. [SJc) . 

In Figs.[SJb) and (d) we plot the peaks' location versus 
inverse system size 1/L for xf and xe, respectively. We 
expect (see for instance the related discussion in Ref. l38l ) 
the following expression to hold for the scaling of the 
peak positions for h c (L) with system size L 



h c (L) =h™ + 



A 



(32) 



where is the result for L — > oo. Data fits give the 
following estimates: = 3.0442(4) and v = 0.625(7) 
[XF, Fig. Mb)} and hf = 3.0442(7) and v = 0.63(1) 
[xe, Fig. EJd)]. We remark that our estimates for the 
location of the QCP are in very good agreement with 
the result from Ref. |3l] and, although quality is lesser in 
this case, our results for v are consistent with the value 
v = 0.6301(8) found in Ref. 

Finally, from the hnitc size scaling analysis performed 
in Sec. Ill CI we expect the following relation to describe 
the behavior of xf on hnite systems in the neighborhood 
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FIG. 5: (Color online) Finite size scaling analysis for the lo- 
cation and height of the peaks in xf [(a) and (b)] and xe [(c) 
and (d)], obtained from the QMC data shown in Fig. [4] In 
panels (a) and (c), the logarithm of the maxima in L" 2 xe and 
L~ 2 Xe, respectively, are plotted as function of InL. Linear 
regression (lines) is applied to the three rightmost data-points 
in each case, yielding the estimates (a) v — 0.623(8) and (c) 
v = 0.615(1) for the correlation length's critical exponent. In 
(b) and (d), the peaks' location h c (L) for, respectively, L~ 2 xe 
and L _2 xe are plotted against inverse system size 1/L. Fits 
(curves) for these results by using Eq. (|32[) yield the estimates: 
(b) = 3.0442(4) and v = 0.625(7) and (d) = 3.0442(7) 
and v = 0.63(1) (the extrapolated values axe indicated by 
the horizontal dashed lines). See main text for details. 



of the QCP 

L- d XF (h,L)=L$- d f XP (L^\h-h c \) , (33) 



and similarly for xe 

L- d X E(h,L) = Li-^f XE (L^\h-h c \) . (34) 

In the above expressions f XF and / XE arc homogeneous 
functions, a priori unknown. Estimates for critical pa- 
rameters can thus be obtained by plotting L _ ~xf and 
£~~ +2 Xe versus L 1 /"^ — h c \ and adjusting the values 
of h c and v until data collapse is achieved. The so ob- 
tained data collapse plots are displayed in Fig. [SJ from 
which we get the following estimates: h c — 3.0440(15) 
and v = 0.625(3) [xf, Fig.^a)], and h c = 3.044(2) and 
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FIG. 6: (Color online) (a) Data collapse for the QMC results 
for L~ 2 \f (a) and L~ 2 \e (b) for the TIM on the square lattice 
and indicated system sizes. Data collapse is achieved for: (a) 
h c = 3.0440(15) and v = 0.625(3) and (b) h c = 3.044(3) and 
v = 0.61(1). In panel (b), data for the smallest system sizes 
are discarded (see main text). 



v = 0.61(1) [xe, Fig. EJb)]. All of these values are in 
agreement with published results \h c = 3.04438(2) from 
Ref. Hi] and v = 0.6301(8), Ref. H?}, but we remark the 
slightly lesser quality of the data collapse achieved for xe- 
Indeed, data obtained from the smallest system sizes in 
Fig. IHb) fail to collapse onto the curve for the largest L 
in Fig. IHl^b) and have not been taken into account when 
performing the analysis. Again, we believe that this is 
explained by the weakly divergent behavior of xe- 



V. DISCUSSION AND CONCLUSIONS 

In summary, we have investigated the scaling proper- 
ties of the fidelity susceptibility xf in the quantum crit- 
ical regime. Large scale quantum Monte Carlo simula- 
tions for the transverse-field Ising model on the square 
lattice, performed by using the scheme introduced in 
Ref. d, confirm the validity of the derived scaling rela- 
tions. Additionally, we also investigate the scaling be- 
havior of the ground-state energy's second derivative, 
d 2 E (g)/dg 2 — — xe(5), a quantity closely related to xf- 



We would like to highlight the fact that the novel QMC 
scheme for computing xf presented in Ref. [1 and dis- 
cussed in detail in the present work opens several re- 
search possibilities within the so-called fidelity approach 
to quantum critical phenomena. Indeed, investigations in 
this field have been so far, to a large extent, restricted to 
one-dimensional systems, while our QMC scheme allows 
for the study of xf for a large class of sign-problem-frcc 
models in arbitrary dimensions. Furthermore, we stress 
that the required modifications in a pre-existing SSE code 
are minimal and that, even though SSE is particularly 
well suited for the task, it is also likely that measure- 
ments for xf(<?) can be implemented within other QMC 
flavors, such as the loop algorithm^ Another potentially 
interesting possibility opened by the QMC method con- 
sidered here involves the study of the finite-T properties 
of xf (which scales as the Bures metric, as shown in 
Sec. llTB) . along the lines of Ref. [H 

A second point worth to emphasize is the particularly 
simple scaling relations for xf derived in Sec. Ill C[ ex- 
pressed solely in terms of the correlation length's critical 
exponent v and that considerably extends the result ob- 
tained by Campos Venuti and Zanardi^ Perhaps even 
more importantly, our scaling analysis does not rely on 
novel concepts such as "quantum adiabatic dimension" 
recently advocated by Gu and coworkers. 40 ' 41 Also, and 
to the best of our knowledge, our Eq. (|2"Tj) is consistent 
with several results for the scaling behavior of xf close 
to second order QCPs presented in the literature, includ- 
ing those compiled in Table I of the review Ref. [! (note 
that some results from the litteraturc quoted in this table 
mistake v for 

Additionally, we also obtain a scaling relation for xe, 
something that allows us to address the important point 
of which of the two quantities, xf or xe, is better suited 
in detecting quantum phase transitions. This question 
is particularly important from the perspective opened by 
the QMC SSE scheme: indeed, as we discuss in detail in 
Scc. lIIII the computational cost for calculating xe can be 
orders of magnitude smaller than the one required in ob- 
taining xf, something in favor of the former as a better 
indicator of quantum criticality, from a practical perspec- 
tive. However, our scaling analysis shows (Sec. Ill Cj) that 
Xe exhibits a weaker divergence (by a factor of z in the 
exponent) than xf, meaning that it might be necessary to 
take into account non-divergent sub-leading corrections 
when performing finite size scaling analysis for xe- Fur- 
thermore, there may even be situations where only xf 
diverges: according to the scaling theory, this happens 
whenever 2/(d + z) < v < 2/d. As a concrete example, 
we mention the case of the QCPs for the CaVO system 
analyzed in Ref. [1, that are preferably detected as a di- 
vergence in xf rather than as a cusp in xe- 

That is, the question of which among xf and xe is 
best tailored to detect an unknown QCP depends on its 
(possibly unknown) universality class and on practical 
matters such as the system sizes that can be reached 
within SSE QMC. From a practical point of view, a pos- 
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siblc strategy consists in evaluating xe on system of up 
to intermediate sizes (where simulations are not too de- 
manding) and search for the presence of peaks hinting at 
a singularity or a cusp in the thermodynamic limit. In the 
affirmative case, simulations for larger systems sizes may 
be performed in order to confirm the occurrence of singu- 
lar behavior for this quantity. If this is not the case, one 
should measure xf(<?) for intermediate sizes and check 
on whether a singularity is more apparent. 

More specifically, we compare now our results for the 
TIM on the square lattice to those obtained, through 
means of exact diagonalizations, by Yu and coworkers. 3 - 
Yu et al. have been able to study xf and xe by con- 
sidering clusters comprising up to 20 sites and have ar- 
rived to the following estimates of critical parameters: 
h c = 2.95(1) and v ~ 1.40. On the other hand, by re- 
sorting on the SSE QMC method discussed in Sec. IIIII 
and on the scaling relations in Sec. Ill C| we are able to 
compute \f for systems with up to N = 28 x 28 sites 
and xe for systems with up to A^ = 48 x 48, arriving at 
the estimates (from the data collapse for xf performed 
in Sec. HVCl) : h c = 3.0440(15) and v = 0.625(3). Our 
estimate for the location of the QCP clearly compares 
much better with results from conventional approaches 
[h c = 3.04438(2) from Ref. HJ than the one found in 
Ref. |32| . And, even more importantly, while our result 
for v is in good agreement with the known result for the 
universality class of the classical Ising model in d = 3 
[v = 0.6301(8), Ref.EzI, the value for v quoted in Ref.|s2 
considerably deviates from it. Again, we suspect that the 
value for v is incorrectly presented as the value for \jv. 

The fact that the analysis employed in Ref. [32| fails 
to obtain critical parameters in agreement with the ones 
from conventional approaches highlights the importance 
of the two main results presented here. First, our SSE 
QMC allows for the computation of xf and xe for much 
larger systems than possible within exact diagonaliza- 
tions, enormously improving the quality of finite-size 
scaling analysis (we remark that results for clusters com- 
prising less than N = 8 x 8 sites are not even taken 
into account in the data collapse performed in Sec. lIVC]) . 
Second, the scaling relations derived in Sec. Ill CI extends 
previous results^ and expresses the scaling dimensions 
for both xf and xe in terms of the correlation length 
exponent. This has the advantage that the exponents 
obtained for xf and Xe can be directly compared to es- 
tablished results for a given universality class, allowing 
us to decide on the validity of the approach. 

We also remark that the scaling relations derived here 
are in agreement with the ones recently derived in the 
field of quantum quenches*^— In this context, the fi- 
delity (and its susceptibility) governs the probability for 
the system to transit to an excited state after a sudden 
change of the coupling constant g away from the critical 
point g c . This expands the range of applicability of the 
concept of fidelity susceptibility beyond the fidelity ap- 
proach to quantum phase transitions^ We might there- 
fore expect that the QMC method presented here, or an 



adaptation thereof, is also applicable in this context. 
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Appendix A: Analytical approximation of an integral 

In this section, we derive useful approximate analytical 
expressions for Eq. (f28|) . 



A(m, n) 



(n-1)! 



m\(n — m — 2)! 



1/2 



First, we note that A(m, n) can be written in a more 
symmetric form 

fn ~\~ 1 

A(m, n) = f(m + 1, n — m — 2), 



where 



f(p,q) 



(p+q + iy- f 1/2 



p\q\ 



/ T P (l-T) 9 dT, 
JO 



We now concentrate on finding efficient analytical ap- 
proximations for f(p, q), in the limit where r = p + q+ 1 
is large. Ultimately, we can apply these estimates to our 
practical case, corresponding to p = m+1, q = n — m—2, 
r = n. 

f(p,q) can be alternatively written 



ji/\ P(1 _ T)qdT 
/ V(l-r)9dT ' 



f(p,q) 



After the change of variable r = | ^1 — ^ p * +g ^ i we obtain 



f(p,q) 



1 - 



p+q 



1 + 



g - v 
I 2 dt 



j -Vp+<i V 



p+q 



1 + 



9-P 

dt 



As a first approximation in the limit r — > 00, keeping 
fixed the ratio 



Xi(p,q) = 



q-p 



(Al) 



we obtain 



f(p,q) 



J o °° e -V+*i(p.«)t dt 
f°° e -T+*i(M)'dt' 

J —00 
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to perform this expansion, we rewrite f(p, q) in the form 



-0.4 



-0.2 
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FIG. 7: (Color online) Results for (a) X(p,q) and (b) f(p,q) 
as a function of p for p + q — 50. In both panels, the thick- full 
line corresponds to the "exact numerical" result, while thin- 
full and dashed curves are respectively first and fifth order 
approximations. The difference between these latter results 
and the exact ones for f(p,q), almost indistinguishable in (b), 
are shown in (c). Notice that the second term in Eq. (|A4[) 
vanishes at p — q [implying f(p, q) = 1/2 for p = q] and the 
approximations become exact at this point. 



which simplifies into 

f(p,q)~G[X 1 (p,q)], 

where 

1 



(A2) 



G{x) 



y 2 /2 



V2^ J- 
1 



dy, 



In order to obtain more systematic and more rigor- 
ous approximate analytic expressions for f(p, q), we quite 
generally introduce A(p, q) such that 



f(p,q) = G[X(p,q)}. 



(A3) 



Since for p — > +oo (at fixed q), one has f(j>,q) — > 0, 
we see that X(p,q) —> — oo in this limit. Moreover, the 
symmetry f(p,q) = 1 — f{q,p) implies that X(p,q) = 
-X(q,p). 

We now look for a systematic expansion of X(p, q) in 
powers of (q — p) and expand the corresponding coeffi- 
cients in (non necessarily integer) powers of 1/r. In order 



f(p,q) 



^(l-^sinh 




dt 


£ (l-t 2 )^ cosh 




dt 



(A4) 

This expression is formally expanded in powers of (q — 
p) and the corresponding coefficients are evaluated for 
large r. This expansion is then matched with the one 
obtained from a similar formal expansion of X(p,q) in 
Eq. (|A3|) . This calculation can be carried out with the 
help of Mathematica, and we finally obtain the fifth 
order expansion in (q — p) of X{p 1 q), which generalizes 
the first order result of Eq. (|A2[) . This expansion can 
be nicely expressed as an expansion in odd powers of 
Xi(p, q), with coefficient having an expansion in integer 
powers of 1/r: 

X{p,q) = a 1 (r)A 1 (p, g ) + a 3 (r)X 1 3 (p,g) + 

a 5 (r)X!(p,q) + -- - , (A5) 

where X\(p, q) is given by Eq. (|A1[) . and with the coeffi- 
cients 

1 19 155 



oi(r) = 1 



12 r 160 r 2 2688 r 3 
1 7 48929 



12 r 360 r 2 362880 r 3 
43 3253 



1440 r 2 362880 r 3 



which were obtained up to third order in 1/r. As sug- 
gested by the above result, one can indeed show that the 
expansion of ct2;+i(r) in powers of 1/r starts at order I. 
Hence, we find that the expansion of Eq. (|A5|) is valid for 
\q — p\ <C r instead of the naive estimate \q — p\ <C ^fr 
which could have been guessed from the quick first or- 
der calculation presented above Eq. (|A2|) . Since one has 
X(p,q) ~ y/r for \q — p\ ~ r, f(p,q) is thus extremely 
close to (p > q) or 1 (p < q) in this regime, with an 
error exponentially small in r. Hence, for all practical 
numerical purpose, it is certainly not a serious problem 
to have an expansion of X(j>, q) limited to \q — p\ <C r. 

We now briefly illustrate the precision of the above 
approximate forms for f(p,q) using the simplest approx- 
imation for X(p, q), given in Eq. (| Al|) . or the fifth order 
calculation of Eq. (|A5[) . For the simplest first-order ex- 
pression of Eq. (|A2[) . the maximal error is less than 10~ 3 
for r > 30, and less than 10~ 4 for r > 275. For the fifth 
order approximation, we find a maximal absolute error 
which is less than 10 -5 , for r > 25, and less than 10~ 7 
for r > 120. 

In Fig. [7ta), we plot X(p,q) as a function of p for 
p + q = 50, for the first and fifth order approximations, 
and for the "exact numerical" X(p, q) obtained after eval- 
uating numerically the defining integral of f(p,q), and 
inverting the relation of Eq. (| A3|) . The plot of the three 
corresponding f(p,q) is presented in Fig. [TJb). 

The maximal error [~ 10~ 3 , see Fig. [7fc)] due to the 
analytical approximations presented in this Appendix is 
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well below our Monte Carlo statistical error and for all 
practical purposes the first-order expression is sufficient. 
Indeed, only in the case of the results presented in Fig. [3j 
we faced the case of r = n < 30. Such small values of 
the SSE expansion order, for which the analytical ap- 



proximations may become not accurate enough, are only 
encountered in the case of very small lattices at high tem- 
perature. In these cases, a simple pre-computation with 
a numerical integration of Eq. (J2SJ) for all values of (m, n) 
can be performed prior to simulations. 
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